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Abstract 



Applying a numerical transfer-matrix formalism, we obtain complex- valued constrained free 
energies for the two-dimensional square-lattice nearest-neighbor Ising ferromagnet below its crit- 
ical temperature and in an external magnetic field. In particular, we study the imaginary part of 
the constrained free-energy branch that corresponds to the metastable phase. Although droplets 
are not introduced explicitly, the metastable free energy is obtained in excellent agreement with 
field-theoretical droplet-model predictions. The finite-size scaling properties are different in the 
weak-field and intermediate-field regimes, and we identify the corresponding different critical- 
droplet shapes. For intermediate fields, we show that the surface free energy of the critical 
' droplet is given by a Wulff construction with the equilibrium surface tension. We also find a 

^ , prefactor exponent in complete agreement with the field-theoretical droplet model. Our results 

extend the region of validity for known results of this field-theoretical droplet model, and they 
indicate that this transfer-matrix approach provides a nonperturbative numerical continuation 
of the equilibrium free energy into the metastable phase. 



1 Introduction 



X. 

^ , Metastability is commonly observed in a wide variety of systems, ranging from supercooled fluids 

and vapors [||, |^ 0] to the electroweak |^, and QCD confinement /deconfinement ^ phase 
transitions, and its description in terms of statistical mechanics has received considerable attention. 
(See e.g. Refs. @, |l|, |ll|, ||, ||] for reviews.) Nevertheless, a fully satisfactory microscopic descrip- 



tion remains elusive |14]. One characteristic feature of metastable phases is that although they 
do not fully minimize the free energy, they nevertheless display long-term stability against small 
perturbations, with lifetimes that can be many orders of magnitude longer than other characteristic 
timescales of the system p5| , p^ , [T7[ |. Based on this observation, numerous efforts to formulate a 
statistical theory of metastability have treated metastable phases as "quasi-equilibrium" phases 
derived from a constrained partition function that excludes or severely reduces the probabilities of 
those microstates that are more probable in equilibrium. In the context of a field-theoretical droplet 
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model with Fokker-Planck dynamics, these ideas led to the result jl^, |T^, ^ that the nucleation 
rate of a metastable phase is proportional to the imaginary part of the analytic continuation of the 
equilibrium free energy into the metastable phase. However, the original derivation of this result is 
only valid for large systems in the limit of ultraweak magnetic fields |19, |2^, and despite extensive 
subsequent studies 



14 21 



, 24, 25, 26, 27, |2^, its domain of validity still remains unclear. 
In one of these subsequent studies |21], it was suggested that the analytic continuation of the 
Ising model free energy could be found from certain analytic properties of the eigenvalues of the 
transfer matrix. Further transfer-matrix work performed on the two-dimensional nearest-neighbor 



Ising model |22, 24] supported this result. Furthermore, these studies found indications for 



the existence of an essential singularity at zero magnetic field, as expected from field-theoretical 
droplet-model calculations |l^, 20, 29, 30, 31, 32]. Since then, this result has been confirmed by 

Monte Carlo simulations 



series expansions 



36, 37], exact diagonalization studies 



4C], rigorous studies using low-temperature Peierls contours ]41], and calculations of two-point 



correlation functions in a bubble model ]M2, 43, 44]. 



Motivated by the transfer- matrix results mentioned above, one of us J^, 46] has introduced 
a "constrained-transfer-matrix" (CTM) method. It generalizes conventional transfer-matrix (TM) 
techniques, providing a means to compute complex- valued "constrained free energies" from the 
eigenvalues and eigenvectors of the transfer matrix. Here we report in detail on the application of 
the CTM method to a short-range-force model, the two-dimensional ferromagnetic nearest-neighbor 
Ising model on a square lattice. We find that the imaginary part of the constrained free energy that 
corresponds to the metastable phase, as obtained by the CTM method, is in excellent agreement 
with predictions of the field-theoretical droplet model. A brief account of some aspects of this work 
has been published elsewhere ]47]. 

One alternative approach to the study of metastability is to perform Monte Carlo simulations 
and measure the lifetime of the metastable phase directly. For the two-dimensional nearest-neighbor 
Ising ferromagnet this was done in Refs. 
|5£ 



4|, H, |5^, 0, |5^, Q and more recently in Refs. ]^, 
|60t] . We compare this lifetime with the imaginary part of the metastable free 
energy as obtained by the CTM method and find very similar dependences on the magnetic field, 
the temperature, and the finite size of the system. This strongly supports the validity of the 
aforementioned proportionality relation between the nucleation rate and the imaginary part of the 
metastable free energy. 

A comparison of the CTM formalism applied to a short-range-force model, as presented here, 
and to two- and three-state models with weak long-range forces, presented elsewhere ] ^5| , |6l| , 
62, 53], clearly reveals the differences and similarities between these models [^^. The analytic 
continuation of the equilibrium free energy in long-range-force models has a branch-point singularity 
at a well-defined non-zero spinodal field ]|l4|, O, 16, 17, I 



55, 56 



58 



27, 28, 45 



61 



|, |6^, whereas 

it has an essential singularity on the coexistence line at zero magnetic field in short-range-force 
models 



29^, 3C 



38, 39 



40, 41, 42, 43, 44]. Moreover, whereas long-range- force 
models exhibit infinitely long-lived metastable phases in the limit of infinite interaction range, 
short-range-force models display finite, albeit very long lifetimes, even in the thermodynamic limit 

|59|, |6C|, |6l, |66l |67|, |68|, |69|. However, the results for 



]22, 48, 49, 5C, 51, 



55^, 



57, 



the imaginary part of the metastable free energy in all of these cases can be explained by assuming 
that the rate-determining step in the decay process is the creation of a critical excitation whose 
free-energy cost can be considered as an activation energy. 

Our results lead to the conclusion that the CTM formalism provides a method to numerically 
continue the equilibrium free energy into the metastable region. In contrast to the analytic con- 
tinuation in field-theoretical droplet models, the CTM formalism does not rely on the explicit 
introduction of droplets. 
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The remamder of this article is organized as follows. In Sec. ^ we explain the CTM method 
and describe the numerical methods used to calculate the eigenvalues and eigenvectors of the TM. 
Section |^ gives a brief overview of the droplet theory with emphasis on those aspects that are 
relevant to the interpretation of our results. The numerical results and their interpretation are 
presented in Sec. ^, followed by a summary and conclusions in Sec. ^. 



2 The Constrained-Transfer-Matrix Method 
2.1 Theory 

We consider the two-dimensional nearest-neighbor Ising ferromagnet on an NxL lattice at a tem- 
perature below its critical temperature T^. The Hamiltonian is given by 

n = -JY.S^Sj-HJ2si, (1) 

where Sj=ibl denotes the spin at site i, J>0 is the nearest-neighbor coupling, the sums run over 
all nearest-neighbor pairs and all sites, respectively, and H is the externally applied magnetic field. 
Throughout this work we use periodic boundary conditions. 

The conventional TM method is a technique to calculate standard thermodynamic state func- 
tions in equilibrium systems, such as the magnetization, the internal energy, the free energy, and 
correlation functions [^0|. In order to define the TM, one first divides the NxL lattice into L 
one-dimensional subsystems of length N. A subsystem configuration \xi) can be written as a direct 
product of single-spin configurations in that subsystem: |x;)=|si^/) • • • |s7v,/)- One then decomposes 
the total Hamiltonian (e.g. Eq. (|^)) into a sum of subsystem Hamiltonians Ti: 

L 

n = Y.n{xuxi+i) . (2) 

1=1 

The elements of the 2^ x 2^ transfer matrix (TM) Tq are then given by 

{xi\Tq\xi+i) = exp[-/57i:(xi,a;;+i)] . (3) 

Here (3=1 /T is the inverse temperature with Boltzmann's constant set equal to unity. Using 
Eqs. (|l|), (^), and periodic boundary conditions in the L direction, the equilibrium partition function 
can be written as 

Z = Tr(T^) =5]A^ (4) 



where Aq, < a < 2 — 1, denote the eigenvalues of the TM in order of decreasing magnitude. Since 
To is finite with all positive elements, the Perron-Frobenius theorem states that the eigenvalue of 
largest norm, Aq, is positive and nondegenerate [^]. Therefore, in the limit L^oo, the equilibrium 
partition function is given by 

Z = A^ (5) 
For the equilibrium free-energy density, /o, we then obtain 

/o = -^lnAo. (6) 

The equilibrium joint and marginal probability densities are given by [70| ] 

Po{xi,xi+k) = {0\xi){xi\{Xq ^Top\xi+k){xi+k\0) (7) 
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and 

Poixi) = {0\xi){xi\0) . (8) 
In terms of the joint and marginal probability densities, the equilibrium entropy is given by 



1 

'n 



^1,^1+1 



Poixi,xi. 



Poixi) 



(9) 



The CTM formalism |45, 46 1 extends this scheme to nonequilibrium states by defining con- 
strained probability densities using the orthonormal eigenvectors \a) that correspond to the sub- 
dominant eigenvalues of Tq: 



Pa{xi,Xl+k) 
Pa{xi) 



{a\xi) {xi I (A„ ^Ta)'''' |xi+fc) {xi+k\a) 
{a\xi){xi\a) . 



(10) 
(11) 



The marginal probability densities, Eq. (11), are normalized and can be interpreted as probability 



densities of subsystem configurations in a constrained state j21, p3, 23, 24, pq|. It is apparent from 



Eqs. (0), d^), ([lO|) , and (11) that the constrained probability densities reduce to their equilibrium 
counterparts for a=0. 

For simplicity, we restrict ourselves in the following discussion to symmetric Tq and assume 
that all its eigenvalues are nondegenerate. (In a large number of cases of physical interest, these 
conditions are or can be made to be satisfied.) The matrices are chosen to commute with Tq 
and can hence be expanded in the eigenvectors, |/?), of Tq: 



Ei/3)/^/3(«)(/5i- 



(12) 



The "reweighted" eigenvalues /i/3(a) are chosen so that T^, for a>0, has Aq, as its dominant 
eigenvalue, rather than Aq- Furthermore, in order for the entire multi-layer system to be in a 
uniform state consistent with Pa{xi), the joint probability densities Pa{xi, xi^k) have to fulfill the 
following requirements |45|: 

1. Stochastic independence at large separation: 



lim Pa{xi,xi+k) 

|fc|— >oo 



Paixi)Paixi+k) 



(13) 



which ensures that fluctuations corresponding to eigenstates orthogonal to \a) decay on a 
finite length scale. This condition is satisfied if and only if |//^(a)|<|AQ,| for all f3^a. 



2. Standard relations between the joint and marginal probability densities: 

E Pa{xuXi^k) = Pa{xi) , 

which is satisfied if and only if iJLp{a)=\a for (5=a. 

3. Non-ambiguity of Eqs. (^) and ( [lT|) for A;=0: 

Pa{xux'i) = Pa{xi)6^^y^ . 



(14) 



(15) 



This condition requires the matrices Tq to be of the same rank as Tq. It is satisfied if and 
only if, for all a, /x/3(a)/0 for all /3 for which Xp^O. 
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These requirements, however, do not uniquely give the detailed form of T^. Since all the A^t^O for 
the system studied here, a reweighting scheme that satisfies Eqs. (|l^), (p^, and (|T5|), and which 
was also chosen in Refs. [H5|, H^, E^, |6l|, E2|, is given by 




(16) 



Note that the Tq, for a>0 are in general not positive matrices. Since only the eigenvector |0), 
corresponding to the largest eigenvalue, can be chosen to have all non-negative elements |^0|, the 
other eigenvectors, which are orthogonal to |0), necessarily contain some negative and some positive 
elements. Therefore, at least some of the matrices Tq, for a>0 may not be positive. 

In order to define "constrained" state functions, it is convenient to decompose the subsystem 
Hamiltonian, 'H{xi,xi^i), into two parts: Tii, containing only interaction terms, and —HM, con- 
taining only terms proportional to the magnetic field. A constrained internal energy per site, Ua, 
and a constrained field energy, —HMa, can then be defined by replacing the equilibrium probabil- 
ity densities in the expression for the equilibrium expectation values of these quantities with the 
constrained probability densities. One obtains 



1 

N 



Pa{xi,Xi+i)ni{xi,Xi+l) 



(17) 



^1,^1+1 



and 



-HMa = -— J2 Pa{xi,xi+i)M{xi,xi+i) . (18) 

^1,^1+1 

Constrained entropy densities Sa are defined in analogy with the source entropy of a stationary, 
ergodic Markov information source (see e.g. JtH): 



1 

'n 



Pa{xi,Xi+i)Ln 



^1,^1+1 



Pa{xi,Xl+i] 



Pa{xi) 



(19) 



Pa{xi,xi+i)Ln{xi\X^^Ta\xi+i) 



^1,^1+1 

Since the elements of T^ are in general not non-negative, this gives rise to an imaginary part in Sa 
through the principal branch of the complex logarithm Lnz=ln |z|+i(/) with z=|z|e"'^ and — 7r<(^<7r. 
Notice also that for a=0, Ua, Ma, and Sa reduce to the corresponding equilibrium quantities. 

In analogy with equilibrium thermodynamics, constrained free-energy densities, fa, are then 
defined by 



fa = Ua-HMa-p-^Sa 



which may also be written as 



-^ln|AQ| + -i- Y Pa{xi,xi+i)Ln 

P P Xl,Xl+l 



(a;;|Ta|x;+i) 

{xi\To\xi 



+1) 



(20) 



(21) 



The first term is analogous to the equilibrium case, whereas the second term vanishes in equilibrium. 
Moreover, the second term is in general complex- valued, since the are not necessarily positive 
matrices. Formally, this second term can be considered as a complex generalization of the Kullback 
discrimination function (see e.g. [^) for Pa{xi,xi+i) with respect to the divergent 'probability 
density' obtained by substituting Tq for in Eq. (|T( 
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2.2 Numerical methods 



In order to calculate any of the constrained free-energy densities in Eq. (pO|), one needs to know 
all of the 2^ eigenvalues and eigenvectors of the equilibrium transfer matrix Tq. The attainable 
system sizes can be increased considerably by making use of the symmetries of Tq (see e.g. [73|). 
Let X be the 2^-dimensional vector space with the basis {j = 1,2, ... ,2^) of subsystem 

configurations. The symmetries of Tq are represented by unitary operators JJ : X ^ X that 
commute with Tq. Applying U to a particular configuration \xj) amounts to a permutation of the 
spins in \xj). The application of U to the configurations in X partitions X into subsets, each of 
which is invariant under U. Since the number N of spins within a subsystem is finite, U has a 
finite period M<N with \J^=1, where I is the identity in X. This means that the eigenvalues of 
U are given by the M roots of unity. Each eigenvalue of U corresponds to one of the subsets of X. 

One can now block-diagonalize Tq by writing it in the eigenvector basis of U, with each block 
corresponding to a different eigenvalue of U. The blocks corresponding to the eigenvalues 1 and 
— 1 are called symmetric and antisymmetric, respectively. They are real and symmetric, whereas 
all the other blocks are Hermitian. Instead of diagonalizing the full 2^ x 2^ transfer matrix all 
at once, one now diagonalizes each block separately, with the submatrices corresponding to the 
Hermitian blocks first converted into real symmetric matrices |74|. 

In our case, Tq is symmetric under cyclic permutations P and reflections R of the subsystem 
configurations. Since P has a period that in general is larger than the period of R, we block- 
diagonalize To with respect to P. In order to identify those states that contribute to the metastable 
phase, we reduce the symmetric block of Tq further by block-diagonalizing it again with respect to 
R. This exploits the fact that the metastable phase has the same symmetry as the stable phase, 
which is symmetric under both P and R |2^, p^ j. 

The diagonalization is performed in two steps. First, we reduce the matrices to tridiagonal 
form by the Householder method, and then we use a QL algorithm with implicit shifts to find the 
eigenvectors and eigenvalues of the resulting tridiagonal matrices [^]. 

Finally, in order to calculate the constrained free-energy densities we transform the eigenvectors 
of the transfer matrix, which are now written in terms of the basis in which Tq is block-diagonal, 
back into the basis {j = 1, . . . 2^) of subsystem configurations. The unitary transformation 

to accomplish this is constructed from the eigenvectors of P and R [ 75 1 . 

It turns out that the imaginary parts of the metastable constrained free-energy densities are 
extremely small (as will be seen later in Fig. 3(b)), although the magnitude of each individual 
term in the sum of Eq. ( ^T]) may be on the order of unity. In order to attain sufficient preci- 
sion, the diagonalization was therefore performed in 128-bit precision on a Cray Y-MP/432 vector 
supercomputer. The total computer time spent in this study was on the order of 1000 CPU hours. 

For temperatures T/J<0.4, the limit on the attainable system sizes is given by numerical 
underflow for small \H\, preventing us from studying system sizes larger than N=9 for T/J=0.4. 
For all higher temperatures studied, the maximum system size attainable is A^=10, limited by the 
available computer memory. 



3 Droplet Theory 

Suppose we prepare the system with H<0, so that its average magnetization is close to —1, and then 
quench the field through H=0 to a value H>0, such that the average magnetization is approximately 
unchanged. This state is then no longer a global minimum of the free energy. Nevertheless, it might 
persist for a period of time many orders of magnitude longer than other characteristic timescales 
of the system |l^, |l^, ^] . This can be used as an operational definition of metastability. 
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Starting with the Becker-Doring droplet model jTo], |Tl|, one assumes that the decay of meta- 
stable phases proceeds through the spontaneous formation of droplets of up spins in the background 
of down spins. Assuming further that the distance between droplets is much larger than their char- 
acteristic size (which is reasonable for sufficiently low temperatures and weak fields fl7\), 
one can treat them as a gas of noninteracting droplets. (In the following, we always assume that 
H is sufficiently weak, so that interactions between droplets are negligible.) The probability of 
creating a droplet is then proportional to the Boltzmann factor e~^^° , where is the free energy 
of the droplet. The final assumption of the Becker-Doring droplet model is that Fj) is given by the 
sum of the surface free energy and the bulk free energy of the droplet. If the surface free energy 
does not change with H and is an increasing and convex function of the volume V of the droplet, 
and if the bulk free energy is a decreasing function of H and V then a critical droplet with volume 
Vc exists so that droplets with V>Vc tend to grow, leading to the decay of the metastable phase. It 
also follows that Vc increases with decreasing H. The quantity Vc can be obtained by maximizing 
Fb- 



3.1 Infinite systems 



We assume for the moment that the system is infinitely large. We then express the surface free 

d-l ^ 

energy, S, in the form 'S=V d E, where d is the spatial dimension, and we denote the difference 
between the bulk free-energy densities of the metastable and the stable phase by A/. Then Fd can 
be written as 

Fb = V^T.-VAf . (22) 



Maximizing Fb yields 



Vr 



(d- 1)S 
dAf 



(23) 



Let Am=mnis — "1st; where mms and rn-st denote the magnetizations of the metastable and stable 
phases, respectively. The approximation Af^\H\Am incorporates the effects of droplet nesting. 



which becomes increasingly important as the temperature increases from zero to Tc |37|. Inserting 
\H\Am and Eq. ( p3[ ) into Eq. (22), the free-energy cost of the critical droplet, Fc{T,H), is obtained 
as 

[d \\H\Am) 



FciT,H) 



(24) 



Within the context of the droplet model, the nucleation rate of critical droplets per unit volume 
is proportional to the probability of creating a critical droplet and can hence be written in the form 

0, g, 0, n 



/(T, H) = Jo(T, ^)e-/^^^(^'^) = Jo(T, H) exp 



(31 



d-l 



with 



'tV fd 



d 



Am 



d-l 



(25) 



(26) 



The proportionality factor Iq{T,H) in Eq. ( p5| ) can be derived by going beyond the Becker- 
Doring droplet model to a field-theoretical droplet model with Fokker-Planck dynamics. Within 
this framework it has been shown |18, 19, 20 1 that for infinitely large systems in ultraweak fields, 
the nucleation rate per unit volume is proportional to the imaginary part of a complex-valued 
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"constrained" free-energy density, Im/ms, obtained by analytic continuation of the equilibrium free 
energy into the metastable phase 20 1: 



I{T,H) = ^\lmU 

TT 



(27) 



Here k is a kinetic prefactor which depends on the details of the dynamics, and the imaginary part 

(28) 



of the constrained free-energy density is given by [^, 31 1 

\lmf^,\=B{T)\H\''exp[-pF,{T,H) [l + 0{H^)]) , 
Inserting (T, H) from Eq. (p^ ) into Eq. (|28| ) , one obtains 



|Im/„ 



B{T)\H\''exp 



1 + 0{H' 



(29) 



where H is given by Eq. (|2q). The function B{T) is expected to be nonuniversal, whereas the 
exponent b is universal and related to surface excitations (Goldstone modes) of the critical droplet 
11,1]: 

J (3 - d)d/2 for l<d<5, d/S 
~ 1 -7/3 for d=3 



(30) 



Omitting the surface excitations in the derivation of Eq. ( p9D leads to 1 75 ] 



(31) 



There is strong numerical evidence that 6=1 for d=2 |l33| , |37| , |57| , |79[, in agreement with Eq. ( pOD 
and disagreement with Eq. (31). 

Under the assumption that the critical droplet is sufficiently compact, one can define its "radius" 
Rc as half of its spanning length along one of the primitive lattice vectors. One can then write 



where fi^ incorporates the detailed shape of the droplet. Using Eq. 

(d- 1)S 



this leads to 



(32) 



(33) 



Making the additional assumption that the critical droplet shape can be obtained from a Wulff 
construction with an anisotropic surface tension, one can relate the shape factor Vld to the surface 
tension. Let cr(n) be the anisotropic surface tension of an interface with normal n. In the Wulff 
construction one first draws a polar plot of the surface tension. According to Wulff 's theorem, the 
shape of the droplet is then given by the inner envelope of the geometric construction obtained by 



drawing a line tangential to the polar plot of o"(n) for each direction n, i.e. |8^, |8l| , 

o-(n) 



\R{f) 



mm; 



n ■ r 



(34) 



Here R{r) is the length of the radius vector in the direction f, and A is a scale factor to be determined 
from the actual size of the droplet. 
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If W denotes the volume bounded by R{f), as given by Eq. (34) with A=l, then the total surface 



free energy of the droplet is given by T.=dW'^/'^V^'^-^^/'^ With S=Sy('^-i)/'=', this leads to 

(35) 




Furthermore, in the Wulff construction with A=l, the length of the radius vector of W in any 
direction is the surface tension in that direction. Therefore W and Vc are related by Vc={Rc/o'o)'^W , 
where cq denotes the surface tension in one of the symmetry directions. Combining this with 
Eqs. (l3^) and (^), one obtains 

so that the radius of the critical droplet from Eq. (Esl) is given by the standard relation 



ft = ^^, (37) 
which, using A/~|i?|Am, can also be written as 

|ii I Am 

The size of a critical droplet therefore increases with decreasing \H\, and its shape changes from a 
square at low temperatures to a circle near Tc. 

3.2 Finite systems 

For A^'^~-'^xcxD hyper cylindrical systems of finite transverse extent N, several distinct field regions 
can be distinguished: an intermediate-field region, where Rc<^N; a small-field region, where Rc^LN; 
and H=0, where both phases coexist. (We do not here discuss the strong-field region, where droplet 
interactions must be taken into account.) 

For Rc<^N, finite-size effects are negligible, so the decay rate of the metastable phase does not 



depend on N. All the results derived in Sec. 3.1 therefore apply in the intermediate-field regime 



For Rc^N, on the other hand, finite-size effects have to be taken into account. This weak-field 
region coincides with the region of validity of finite-size scaling at first-order phase transitions 
described in Refs. |8|, pi Isl. 



We restrict the following discussion to d=2. It has recently been shown rigorously 



90, 91| that for very low temperatures and \H\/J<4:, a large metastable Ising ferromagnet in two 
dimensions with local Metropolis dynamics decays through a single nucleating droplet of the shape 
shown in Fig. 1(a). This is an (lc—l)xlc rectangle of overturned spins with a single additional 
overturned spin as a "knob" on one of its long sides. The length lc=\2J/\H\'] is the smallest integer 
larger than 2J/\H\, where 2J/\H\=limT~,o 2Rc- Flipping the remaining lc—1 spins on the side with 
the "knob" results in a net gain of energy, and the resulting IqXIc droplet is the low-temperature, 
discrete-lattice equivalent of the critical droplet discussed in Sec. |3.1| . Finite-size eff^ects due to the 
possibility of subcritical droplets wrapping around the lattice in the transverse direction are to be 
expected for or (extrapolating in an approximate fashion to nonzero T) 2Rc>N—2, which 

yields \H\<2ao/[{N-2)Am]<2J/{N-2)=H2. 

However, critical droplets of equilibrium shape are not the only configurations comparable to 
the system size that may nucleate the stable phase by wrapping around the lattice. It is easy to 
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Figure 1: Sketch of two critical excitations that contribute to the decay of the metastable phase 
at weak positive fields and low temperatures. Configurations shaped like the one depicted in 
(a) dominate for \H\^H2=2J/{N — 2). The "rod"-like configuration shown in (b) prevails for 
\H\<Hi=2J/{N - 1). See details in Sec. |]|. 

show that for \H\<H2, still lower-energy configurations that are capable of spanning the system 
by flipping only one single spin can be constructed by successively removing one layer of + spins 
from one of the short sides of the droplet shown in Fig. 1(a), until one reaches the (iV— l)xl "rod"- 
shaped cluster shown in Fig. 1(b). For H2>\H\>2J/(N—l)=Hi, a further energy reduction can be 
achieved by adding (A^— l)xl slices back onto one of the long sides of the cluster in Fig. 1(b), until 
one reaches an (A^— 1) x (iV— 1) cluster. This does not hold true for \H\<Hi. Here the "rod"-shaped 
cluster of Fig. 1(b) is the lowest-energy configuration that can nucleate the stable phase by flipping 
a single spin. 

The above observations lead to the following predictions for the finite-size scaling of |Im/jns|- 
For \H\>H2, finite-size effects should be negligible, for Hi<\H\<H2 the behavior should be rather 
complicated due to the large number of competing clusters with nearly degenerate energies, and 
for \H\^Hi the behavior should be determined by "rod"-shaped clusters like the one shown in 
Fig. 1(b). 

To leading order, the free-energy cost Fc of creating the "rod" -shaped critical cluster in Fig. 1(b) 
is given by its surface free energy only, so that Fc^2a{T, H)N. In analogy with the equilibrium 
surface tension, a{T, H) here denotes the "surface tension" of the critical droplet at non-zero fields. 
Inserting Fc into Eq. (p8[), one obtains 

|Im/^3| oc e-2/^'^(^'^)^ . (39) 

The largest corrections to Eq. ( |3^ ) should include a bulk term in the exponential, which is of 
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the form —2\H\{N — 1) for the particular critical excitation considered here, and an A^-dependent 
prefactor. 

In the infinite-system limit H2^0- We then expect that a{T,H) -^aeq{T) for large A^, where 
(7eq(T) is the exact equilibrium surface tension as obtained from Onsager's solution of the two- 
dimensional Ising model |]8T| , |85| . 

At H=0 one can show (see Appendix^) that |Im/ms| as defined in the CTM formalism, Eq. (pOD, 
is given by 

|Im/^s| oc {mN)-' . (40) 
Here, ^Ar=ln |Ao/Ai| denotes the correlation length which at H=0 is [|^, 86, 37| 

Cn oc Arl/2g/3..,(T)7V ^ 



S6 , 



Inserting Eq. 



(41) 
(|1]) into 

(42) 



with a nonuniversal proportionality constant 
Eq. (^0|) one obtains 

|Im/^s| oc Ar-3/2g-,3..,(T)7V 

Note that the coefficient of cJeq in Eqs. ( |39[ ) and (42) differ by a factor of two. Physically, this 
corresponds to the fact that, whereas the rod-like droplets that dominate for \H\^Hi correspond 
to pairs of tightly bound transverse interfaces, the interfaces are unbound for H=0 [B6|. 



4 Numerical Results and Discussion 



Since the metastable phase has the same symmetry as the stable phase [53, 23, 24], which is 



symmetric under cyclic permutations and reflections, we consider only those eigenstates of the 
TM that are symmetric under cyclic permutations and reflections. As an example, in Fig. 2 we 
plotted the eigenvalue spectrum as —liiXa/{PN) (Fig. 2(a)), which corresponds to the first term 
in Eq. (|l]), and the constrained magnetizations (Fig. 2(b)) for T/J=1.0 (r/rc=0. 44068 . . .) and 
N=8. When the eigenvalue spectrum is plotted as in Fig. 2(a), the branch corresponding to a=0 
is the equilibrium free-energy density, Eq. (^. In Fig. 2(b), it corresponds to the branch with 
magnetization close to +1.0. 

As seen in Fig. 2(a), for weak magnetic fields, certain eigenvalues group together into fans. In 
order to identify which eigenvalues correspond to which subsystem configurations, we calculated for 
T=0 the internal energies, Ua, and the magnetization, M^, for all those subsystem configurations 
that are symmetric under cyclic permutations and reflections. Configurations with the same number 
of interfaces parallel to the transfer direction have the same value of Ua- We also find that at H=0 
and r>0, the averages calculated over all those — In Aq,/(/5A^) that belong to the same fan are 
nearly equal to the T=0 values of Ua- (For example for A^=10 and T/J=1.0, the discrepancies 
between Ua and the corresponding averages are smaller than 2%.) We therefore conclude that 
all the eigenvalues that belong to the same fan correspond to subsystem configurations with the 
same number of interfaces. The absence of exact degeneracies between the eigenvalues that belong 
to the same fan for T>0 is an entropy effect. Finally, we find that away from the crossings, the 
negative of the slopes of the — In Aq,/ {I3N) are approximately given by the magnetizations Ma (with 
discrepancies that are again smaller than 2%). 

the 



In accordance with earlier transfer- matrix studies |21, 23, 25, 61, 32 



92| 



eigenstate \a) corresponding to the metastable phase at a given field was identified as the one with 
the largest magnetization opposite in direction to the magnetic field (marked by a thick solid line in 
Fig. 2 in the field region where this — In Aq, is not the uppermost branch). As was already observed 
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previously |^2|, such an eigenstate can be identified. Furthermore, the metastable eigenvalues vary 
little with A^, indicating the existence of a well-defined thermodynamic limit for that state |22]. 

At certain values of the field the metastable branch seems to cross different eigenvalue branches. 
These are in fact avoided crossings, which only become exact crossings at T=0. Consequently, 
the metastable branch does not consist of only one eigenstate, but is composed of a succession 
of eigenstates at different fields. For r=0, we can use the relation — In Xa/{PN)=Ua—MaH to 
calculate exactly the value of the magnetic field at which the m*^ crossing involving the metastable 
branch occurs: 

2J 

Ho = 0, and Hm = , with m=l, . . . , iV-1 . (43) 

N — m 

For the non-zero temperatures studied here, the field values of the avoided crossings involving the 
metastable branch deviate by at most 3% from their values at T=0. For m=l and m=2 Eq. (43) 



corresponds to the definitions of Hi and H2 given in Sec. 3.2 



From Eq. (|4^), the number of avoided crossings in any neighborhood of H=0 is seen to go to 
infinity as A^— >oo. This is consistent with the existence of an essential singularity in the free energy 
at H=0 in the thermodynamic hmit |l9|, |2|, |3|, ||, ||, ||, |3^, ||, ||, |4^, ^, |4 1 . 



In particular, it is reminiscent of the results of recent studies which found an essential singularity 



in the susceptibility at H=0 p^ , 43, 44|. In that work the susceptibility was calculated in an 
ensemble in which the magnetization was restricted to negative values for positive fields. It was 
found that in the limit H^O, the singularity manifests itself through an infinite number of poles on 
the positive real axis, whose locations are given by H(xl/n, n=l, 2, . . . , 00. However, there is specific 
disagreement between the values of the fields at which these poles occur and those corresponding 
to the avoided crossings in the TM spectrum. Since the TM eigenvalue crossings coincide with 
singularities in the metastable lifetimes recently found in rigorous low-temperature calculations 



p8| , we believe the numerical discrepancy with Refs. |4^, 43, Q may be due to approximations 
used in that work. 

Figure 3 shows the real parts (Fig. 3(a)) and the imaginary parts (Fig. 3(b)) of the constrained 
free-energy densities for T/ J=1.0 and A^=8. The composite branch corresponding to the metastable 
phase is marked by thick solid lines in both panels. The general features of the stable and metastable 
phases in Fig. 2(a) are reproduced in Fig. 3(a). For the real part of /q, the composite metastable 
branch is, except near the crossings, almost identical to the composite metastable branch in the 
eigenvalue spectrum, shown in Fig. 2(a). In the spectrum of the imaginary parts, the composite 
branch corresponding to the metastable phase consists of different lobes, each of which corresponds 
to a different eigenstate. The crossings for the metastable llm/^l correspond to avoided crossings 
in the eigenvalue spectrum of Fig. 2(a). Split lobes (such as the one at i//J«0.8) occur when 
the metastable branch intersects branches from fans corresponding to more than two interfaces. 
Notice the extreme smallness of the values of the minima of the metastable llm/ol, especially for 
weak magnetic fields, and how their range extends over roughly ten decades. We also observe, in 
accordance with our expectations for the finite-size scaling behavior of |Im/ms| outlined in Sec. 3.2 , 
that the value of the metastable |Im/a| at the first minimum is approximately given by the square 
of its value at H=0. 

In order to avoid complications introduced by the mixing of two or more eigenvectors in the 
vicinity of the crossings (where, for example, the constrained magnetizations of the metastable 
eigenstates deviate appreciably from —1), we concentrate on the minima of llm/^l, which are 
located away from the crossings. We denote the values of the fields at which these minima occur by 
-ff™™, such that Hm-i<H'^'^<Hm for m=l, . . . A^— 1 where Hm are the crossing fields of Eq. (43). 

Figure 4 is a semi-log plot of only the metastable llm/^l vs. inverse field, J/H, for A^=9 and 
A^=10 at r/J=1.0. The thick straight line was drawn through the two minima for iV=10 between 
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Table 1: The diameters of critical droplets, 2Rc{T,H), evaluated at (a) Hf^'iN), (b) Hf^'iN), 
and (c) H^^^i{N). Here N is the largest system size used at a particular field and temperature. 
At H=Hf^^_^^{N) and for T/J>1.1, the minimum of the metastable |Im/Q,| was so shallow that we 
were unable to locate it. 
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J/H=3.0 and J/H=4.0, but, in accordance with Eq. (29), it also connects minima at smaller values 



of J/H for both N=9 and A'^=10. In agreement with our expectations outlined in Sec. 3^ above, 
we find entirely different dependences on N and H for the metastable llm/d in the weak- field and 
the intermediate-field regions. As is apparent from Fig. 4, for if>i?™° llm/al is independent of N. 
For ff<iJ™™, however, |Im/Q,| depends strongly on N. This is further supported by Table [l|, where 
we compare the diameter, 2Rc, of the critical droplet with the system size, A'^, at three different 
values of the field, ff™™, ff™™, and /f™™]^. We find that for all temperatures studied, the size of 
the critical droplet exceeds A^ at ff™™, is comparable to N at i^™™, and is close to unity at H^^^^. 
The quantity Rc was calculated from Eq. ( |38| ) using for ctq the exact equilibrium surface tension, 
crcq(T), as obtained from Onsager's solution of the two-dimensional Ising model ]8l| , 82, 85]. We 
also set Af=2meq{T)H, where 



eq 



(T) = (1 - cosech^2/3J)i/s (44) 



is the exact zero-field equilibrium magnetization ||93[]. These identifications will be justified below. 

The different dependences of the metastable llm/d on A^ and H are discussed in more detail 
in the two subsections below. 
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4.1 The weak-field region 

Figure 5 shows a semi- log plot of the metastable A^llm/ol (plotted as *) and (shown as o) vs. 
N at H=0 and T/J=1.0. The dashed straight lines are guides to the eye drawn through the points 
at A^=9 and A^=10. The fact that these lines also approximately connect points at smaller values 
of N indicates that, in accordance with Eqs. (^Tj) and (|4^), the metastable llmf^l and vanish 
exponentially with N. 

The surface tension was calculated from either Eq. ([4l| ) or Eq. ( p2|) by taking either the difference 
of the logarithms of the metastable |Im/Q,| between two successive N or the difference of the 
logarithms of between two successive A'^. Taking into account appropriate correction terms 
arising from the A^-dependent prefactors led to estimates for the surface tension that agree very 
well with (Teq(T), for both the metastable |Im/Q,| and ^'j^^. Except for the lowest temperature 
studied (which is hampered by numerical underflow), the differences between our estimates and 
aeq{T) are smaller than 2%, for both the metastable |Im/Q,| and 

For the weak- field region. Fig. 6 shows a semi- log plot of the metastable llm/^l at vs. 
for all the temperatures studied. In accordance with Eq. (p9|), we find that the metastable |Im/Q,| 
vanishes exponentially with N. The slopes of the straight lines in Fig. 6, divided by —2/3, give an 
estimate for a(T, H) at the different temperatures. 

We calculated a(T,H) by taking the difference of the logarithms of the metastable |Im/a| in 
Eq. (|3^) at two different system sizes. In order to eliminate noticeable effects due to even-odd or 
odd-even pairs of system sizes, cr(T, H) was obtained by using only odd-odd or even-even pairs of 
A'^. According to Eq. (^), the value of the field. Hi, where the first crossing occurs, converges 
to zero in the limit N^oo. Therefore, i?f™<ffi also converges to zero, so that in the infinite- 
system limit, a(T,H) should converge to the equilibrium surface tension creq(T). Assuming linear 
convergence, three different values of a{T,H), as obtained from Eq. (l39|) , are required in order to 
extrapolate a{T, H) to the infinite-system limit. Again, to avoid effects due to even-odd or odd- 
even combinations of system sizes, only the cr(T, H) calculated from system sizes with the same 
parity were used to extrapolate cr(T, H) to the infinite-system limit. The resulting extrapolations 
were then averaged to obtain an estimate of cr(T) in the infinite-system limit. This procedure could 
not be carried out at T/ J=QA due to numerical underflow for A^>6. Instead, the value of cr(T) 
quoted at this temperature corresponds to (j{T) as calculated directly from Eq. ( |39| ) using A^=4 
and 6. 

In Fig. 7 we compare (t{T) obtained from Eq. ( |39| ) (shown as x) with the exact equilibrium 
surface tension, creq(7') (shown as the solid curve), for all the temperatures studied. Except for 
the lowest temperature, the agreement between the extrapolated result and creq(T) is excellent, 
indicating that our assumptions about the significance of llm/^l and about the shape of the critical 
excitation in the weak- field region are correct. It also suggests that in the thermodynamic limit the 
surface free energy of the critical excitation is given by the exact equilibrium surface free energy. 

We tried to improve the estimates of (j{T) by including various corrections. As an example, 
in Fig. 7 we also plot extrapolated results obtained from incorporating a bulk term of the form 
—2H{N — 1) in the exponential in Eq. ( |39|) (shown as □). The fact that we do not find universal 
improvement over the results obtained without corrections indicates that other terms besides a 
bulk term might be equally important. However, including various A^- or //-dependent prefactors 
into Eq. ( p9| ) (as in Eq. (41) or Eq. (42)) separately or in conjunction with the bulk term did 



not consistently improve our previous results. These findings, together with the small difference 
between all our estimates, imply that the surface term in the free energy of the critical excitation 
is the dominant one. 
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Figure 8(a), where H^^™ is plotted vs. N for various temperatures, shows that i?™™ scales as 

Hf''/J = ci(r)iV-°(^) . (45) 

Here ci(T) is a proportionality factor. The temperature dependence of the exponent a{T) is shown 
in Fig. 8(b). It is approximately two for r/J=1.0 and 1.2 and drops to a value below 1.5 for 
r/J=0.4. Since H™™<Hi^ according to Eq. (^), a(T)>l for all temperatures. We also expect 
\\m.x^Tc <^(^)=2, corresponding to a correlation length that grows linearly with N at the critical 
temperature, in agreement with critical finite-size scaling theory |^3|. However, neither the origin 
of the N dependence of i/™"^, nor the detailed temperature dependence of the exponent Oi{T), is 
apparent to us at this time. 



4.2 The intermediate-field region 



A very different scaling behavior is found for i?™™^i?^2J. In this field interval. Fig. 4 indicates 
that the logarithm of the metastable |Im/Q,| is independent of A^. Only at fields do we find 

some remaining finite-size effects, which are i;3%. At stronger fields they vanish almost completely, 
becoming smaller than 0.1%. If, as we propose, the metastable llm/d is related to the free energy 
of a critical droplet by Eq. (|28|), then this field interval of A-independent llm/d should correspond 
to critical droplet sizes between N and unity. This was already seen in Tables ^(b) and |l|(c). For 
all temperatures studied, we found that the diameter of the average critical droplet ranges from 
about eight lattice units for A=10 near the crossover field i?™™, down to one lattice unit for all A^ 
near H/ J=2. 

It can already be inferred from the thick straight line in Fig. 4 that the leading-order field 
dependence of the metastable |Im/a| at intermediate fields can be described by Eq. (29) with an 
approximately field-independent value for the parameter S. The slope of this line gives a rough 
estimate for H, which was defined in Eq. (|2^ . Below, we emphasize the lack of an H dependence 
, S, (To, and meq by using the explicit notation H(T), S(r), cro(r), and meq(T). 



m 



In order to determine H(r), b, and B(T) in a more systematic way, and thereby to study the 
validity of Eq. (^) for intermediate fields in greater detail, we performed linear least-squares fits 
on the logarithm of the minima of the metastable |Im/Q,|, excluding terms of order higher than 
two in the exponent. That leaves four independent parameters. However, assuming that B{T) 
is independent of the field, it can be eliminated by differentiating ln|Im/ms| with respect to the 
inverse magnetic field. With |Im/ms| as given in Eq. ([29|) and d=2, one obtains 



din |Im/n 

~W/\H\: 



pE{T) + b\H\ + 0{H^ 



(46) 



In order to utilize Eq. (^), we use the two-point finite-difference estimate for the derivative of 
ln|Im/ins| based on the metastable |Im/„| at two successive minima, -fT™™ and -f^^™i: 



A(ln|Im^|) 
A{l/\H\) 



\n\lmUiK 



mm\ 
m ) 



ln|Im/,(i/-V"i: 



ii\m 



l/\H. 



ni+l I 



/?H(r) + hH,s + 0{H^) 



(47) 



with 



cflf 



TTvam trmm 



if" 



In 



tjmm 



(48) 



We then perform three-parameter fits on A(ln[Im/Q])/A(l/|/7|). 
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In the following, all the results reported for H(T) and b were obtained from three-parameter fits, 
whereas the results given for B{T) were calculated from two-parameter fits with H(r) and b fixed. 
It should be emphasized, however, that since our results for llm/^l are numerically exact, we are 
using the least-squares method merely as a tool to determine the above parameters numerically. 
We have therefore assigned error bars such that P^r degree of freedom is approximately unity. 

Typical results of the three-parameter linear least-squares fits are shown in Fig. 9. The ratio 
— A(ln |Im/a|)/A(l/|//|) vs. H^g is plotted at four different temperatures. The symbols represent 
the results obtained with the CTM method, except at zero field, where we plotted the numerically 
exact /3Heq(r)=/3i;2q(r)/(4Am), with Am=2mcq(r). Here SeqlT) was obtained by integrating 
numerically over the equilibrium droplet shape, which in turn was calculated by combining the Wulff 
construction with Onsager's exact zero-field anisotropic surface tension The equilibrium 

magnetization meq(T) is given by Eq. (Q) above [^]. The solid curves represent our fits. The 
dashed straight lines start at the exact /3Heq(T) at zero field and have a slope of unity, corresponding 
to the value of 6=1, expected from field theory, Eq. (30). It should be clear from Fig. 9 that our 
data in conjunction with the exact result at zero field can not be fitted by straight lines which 
would only involve b and H(T) as free parameters. The reason lies in the competition, in Eq. (p9|), 
between the 0{H^) correction term in the exponential, which dominates for strong fields, and the 
prefactor, which becomes important for weaker fields. The accessible fields are sufficiently strong 
that the 0{H^)-ierm is crucial and has to be included into the fits. Furthermore, the field at 
which the domination changes from one term to the other decreases as the temperature is lowered. 
That means that for lower temperatures, the range of fields we can cover does not extend as far 
into the region of weaker fields where the power-law prefactor dominates, as it does for the higher 
temperatures. This makes it increasingly difficult to obtain reliable estimates for the prefactor 
exponent 6 as T is lowered. In the intermediate-temperature region of 0.9<T/J<1.1, however, we 
find excellent agreement between the exact /3Heq(r) at zero field and the corresponding quantity 
obtained from our fits. 

This agreement is further illustrated in Fig. 10. It shows S(T) divided by its equilibrium 
value, Heq(T), vs. T, as obtained from three-parameter linear least-squares fits, including those 
shown in Fig. 9. For all T/ J>0.4 these fits were performed on data corresponding to A^=10. Due 
to numerical underflow, we were limited to a maximum system size of A^=9 and flelds i?>i?™™ 
for r/J=0.4. At temperatures r/J>l.l, on the other hand, the minimum at -ff™!!^! becomes so 
shallow that we were unable to locate it. That reduces the number of degrees of freedom in the 
fits for the higher temperatures, T/J>1.1, by one as compared to the temperatures 0.6<r/J<1.0, 
making the fits for these temperatures slightly less reliable than the others. Two sets of data are 
shown. The data set corresponding to the diamonds was obtained from fits using all the discernible 
minima of the metastable |Im/a| in the interval H2^^^<H<2J. The data set corresponding to the 
crosses with error bars was obtained by excluding the minimum at (For clarity, the error 

bars for the diamonds, which are comparable to those for the crosses, are not shown.) This was 
done in order to illustrate the infiuence of uncontrollable finite-size effects due to competing critical 
excitations present in the metastable llm/^l at ff™™. The largest finite-size effects occur at the 
higher temperatures T/J> 1.1. Corrections due to finite-size effects stay within the limits of the 
error bars for T/J=0.8, 0.9, and 1.0 and increase for both higher and lower temperatures. In the 
following, we always exclude from the fits the minimum of the metastable llm/ol that occurs at 

The fits shown in Fig. 10 indicate that our results for H(r) are in excellent agreement with the 
hypothesis that the surface free energy of critical droplets is equal to that of equilibrium droplets of 
the same size at the same temperature, and that Am equals twice the exact zero-field magnetization 
meq(T). Except for T/J=l.l, we find agreement to within 10% between the estimated H(T) and 
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Heq(r). 

Independent confirmation of our results has been obtained by using Eq. (27) in conjunction 
with direct measurements of the hfetime of the metastable phase in Monte Carlo simulations |57, 



5|, H, |0|. 

If we set H(T) to its equilibrium value Heq(T), thereby reducing the three-parameter linear 
least-squares fit to a two-parameter fit, we obtain estimates for the exponent 6, which are plotted 
in Fig. 11. As in Fig. 10, all the estimates were calculated for H^^^<H<2J and A^=10, except at 
r/J=0.4, where we used A^=9. For T/J>0.8 our results are consistent with b=l. Because of the 
aforementioned competition between the prefactor term, which involves the exponent b, and the 
0{H^) correction term in the exponential of Eq. (p9D, the estimates for h become less reliable at 
lower temperatures. As a result, we do not believe the large deviations of the estimate of h from 
unity for T/J<0.8 to be physically meaningful. The general agreement of our estimates with 6=1, 
in agreement with Eq. ( |30|) and in disagreement with Eq. (|3T|), confirm the importance of surface 
excitations on the critical droplet in determining the prefactor exponent h. Our results are fully 
consistent with previous numerical estimates [^, |5^, ^ . 

Figure 12 shows a comparison of Heq(T), indicated by the solid line, with estimates for H(T) 
obtained from a two-parameter linear least-squares fit, in which b is fixed to 1. The droplet shape 
varies between a square at r=0 and a circle at T=Tc (both shown for the whole temperature 
range by dashed lines). As in Fig. 10, we find close agreement between S(r) for the critical 
droplet with the corresponding value Heq(T) for the equilibrium droplet. Together with Fig. 10, 
these results strongly support the notion that at a given temperature, the metastable phase decays 
through the formation of critical droplets having the same shape as equilibrium droplets at the 
same temperature. This also justifies the use of the exact equilibrium surface tension to estimate 
the size of the critical droplets in Table ||. 

We also determined the prefactor B{T) by reducing the four-parameter linear least-squares 
fit to Eq. (29) to a two-parameter fit by setting H(T) to its equilibrium value and 6=1. Our 
results are displayed in Fig. 13 as x with error bars. They are consistent with an exponential 
dependence of B{T) on l/T, shown by the dashed line which was obtained from a weighted linear 
least-squares fit to our results. This temperature dependence of B{T) can be obtained by adding 
an approximately temperature independent term AF to the free energy of the critical droplet. Our 
linear least-squares fit indicates that AF~1.8. This result is supported by a recent Monte Carlo 
study Using continuum droplet theory, values of B{T) were calculated from measurements 

of the lifetime of the metastable phase and are shown as O with error bars in Fig. 13. They also 
follow a straight line with a value of AF~1.2 which is of the same order of magnitude as the CTM 
result. The discrepancy between the two estimates of AF is likely due to a temperature dependent 
kinetic prefactor in the Monte Carlo results. Also shown in Fig. 13 are values of B{T) as calculated 
from the Becker-Doring droplet model by Harris |3^. Clearly, the Becker-Doring droplet model 
is inadequate to describe the temperature dependence of our results. An additive term in the free 
energy of the droplet was also found by Jacucci et al. |39] which they interpreted as a curvature- 
dependent surface free energy. However, although their additive term is field independent, it does 



depend on the temperature. Recent rigorous studies predict AF=AJ at very low T pH, n9, BO, 91 



We note that our attempts at providing a simple, heuristic extrapolation of this result to the 
temperature range covered by the CTM and Monte Carlo studies by using AF(T)~2(t(T) were not 
successful in fitting the numerical results. A clarification, especially of the physical origin of AF 
as found in the CTM and Monte Carlo results, must await further study. 

For the interpretation of the results of the CTM method in terms of the droplet model to 
be consistent, the critical droplets have to be compact with well defined boundaries. Explicit 
verification of this is given in Appendix |H. 
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5 Summary and Conclusions 



We successfully applied the constrained-transfer-matrix method to the two-dimensional 

nearest-neighbor Ising ferromagnet in an external magnetic field to obtain complex-valued con- 
strained free energies. In very good agreement with nucleation theory, we found several different 
field regimes, set apart by characteristic field and system-size dependences of the metastable |Im/Q,|. 
At H=0, we find that the metastable |Im/Q,| is proportional to {N^]\f)~^ and therefore decays ex- 
ponentially with N. For fields 0<H<Hi, where the size of the critical droplet as determined from 
droplet theory is larger than the transverse system size A^, the decay of the metastable phase pro- 
ceeds through the formation of "rod"-like excitations that extend (A'^— 1) lattice units in the finite 
direction and are only one lattice unit wide. Introducing a "surface tension" for nonzero fields, this 
again leads to an exponential decay of the metastable llm/d with N. Due to the binding of two 
interfaces in the rod-like excitation, the argument of the exponential factor in the metastable |Im/Q,| 
at 0<-ff™™<-ffi is approximately twice as large as the one at H=0. Extrapolating the nonequilib- 
rium surface tension to the infinite system limit, we find very good agreement with Onsager's exact 
equilibrium surface tension. For H>H2, the decay of the metastable phase is dominated by critical 
droplets that are smaller than the finite width of the system. The metastable |Im/Q,| in this field 
regime is therefore A- independent. The functional form of the imaginary part of the metastable 
free energy is the same as the one obtained from field-theoretical droplet models for infinitely large 
systems in ultraweak fields |3^ (see Eq. (p9D). Also, our numerical results indicate that the crit- 
ical droplets have the same surface free energy as equilibrium droplets of the same size at the same 
temperature. This and our estimates for the prefactor exponent b are in excellent agreement with 



series expansion |33, 34, 35, 3£, exact diagonalization [^], and Monte Carlo results [39, 40 1. 



A residual A-dependence remains in the transition region Hi<H<H2, where a number of critical 
excitations of different shapes have comparable energies. Finally, in Appendix |^, we have given 
some estimates for the field /^mfsp beyond which standard droplet theory is expected to become 
unreliable. For iJ>iJMFSP> we could not resolve a metastable imaginary free energy from all the 
other branches of |Im/Q,|. 

Our results for the dependences of the metastable |Im/Q,| on system size, field, and temperature 
can be compared with direct measurements of metastable lifetimes from Monte Carlo studies Hq, 



49, 50, 56, 57, 58, pO|. The picture that emerges from Ref. [^] involves four distinct field 

regions, in which the decay proceeds via different mechanisms. For very weak fields, the lifetimes 
obtained from the Monte Carlo simulations grow exponentially with the area of an interface cutting 
across the system, indicating that the excitations that lead to a decay of the metastable phase are 
of a size comparable to the finite dimension of the system. For intermediate fields, the lifetime 
is proportional to the inverse system volume, indicating that a small, size-independent nucleation 
rate is the rate-determining step. This "single-droplet region" |5^] is characterized by decay via a 
single critical droplet. For yet stronger fields, the lifetime is size independent, indicating that the 
rate-determining process is the nucleation and growth of a finite density of droplets. This region 
was called the "multi-droplet region" in Ref. |57|. Finally, in the "strong- field region" beyond 
-ffMFSPj the droplet picture becomes inadequate. 

In the CTM study presented here, the "weak-field region", where |Im/a| depends exponentially 
on the system size, and the "intermediate-field region" correspond directly to the weak-field and 
the single-droplet regions found in the Monte Carlo study in Ref. |57|. In the CTM study the 
"weak-field region" extends from H=0 to approximately the field of the first eigenvalue crossing at 
Hi /J and the "single droplet region" extends from approximately Hi/ J to approximately -ffMFSP 
given in Appendix with no clearly distinguishable "multi-droplet region" . The values for -ffMFSP 
obtained in the study presented here are of the same order of magnitude as the ones obtained by 



18 



Monte Carlo simulation [p7| . 

Comparing the field and temperature dependence of |Im/a| with that of the metastable lifetime 
from Monte Carlo calculations [57, |5^, pO| shows that for H<Hmfsp both quantities can be 



obtained to leading order from the same droplet theory. The only difference is a kinetic prefactor 
that enters into the calculation of the lifetime pO], which depends on the particular dynamic 



chosen in the Monte Carlo algorithm [57|. This is in complete agreement with and therefore supports 
the validity of Eq. ( |29| ) not only for ultraweak fields but also for intermediate fields. 

Comparison of the results obtained for the typical short-range-force model, presented here, with 



CTM studies of two- and three-state systems with weak long-range interactions |45|, |61|, |6^, |63[ 
clearly shows the differences and similarities in the properties of the metastable phase in the two 
types of systems The lifetime in short-range-force models becomes large but stays finite even 
in the thermodynamic limit, with an imaginary part of the metastable free energy that displays 
a droplet-model type essential singularity at zero magnetic field. Long-range-force models, on 
the other hand, exhibit metastable phases that become infinitely long-lived as the range of the 
interactions goes to infinity fl^, |l6|, with an imaginary part of the metastable free energy that 
exhibits scaling behavior consistent with a well-defined non-zero spinodal field and a corresponding 
branch-point singularity |6^, As the CTM studies for those systems show, the region in which 
the metastable |Im/a| decays exponentially with extends all the way out to the exactly known 
mean-field spinodal. Nevertheless, in both types of models, the leading-order dependence of the 
metastable |Im/a| on the magnetic field and the temperature is, over a wide range of fields, given 
by a Boltzmann factor involving a well-defined free energy of formation for a critical excitation in 
the metastable phase |6^, |63| . 

All our results strongly support the conclusion that the constrained-transfer-matrix formalism 
provides a nonperturbative method to numerically continue the equilibrium free energy into the 
metastable phase that, in contrast to the analytic continuations obtained in field-theoretical droplet 
models, does not rely on the explicit introduction of the particular excitations through which the 
metastable phase decays. 

A question of current interest is whether or not the decay of metastability in anisotropic square- 



lattice Ising ferromagnets [94, ^] proceeds through equilibrium-shape critical droplets [81 



A study of this topic by the CTM method is in progress [96| 
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A The relationship between |Im/ins| and at H=0 

As may be seen from Fig. 3, the metastable phase for \H\<Hi is represented by the eigenstate |1), 



so that Im/ms=Im/i at H=0. Using Eqs. (|l^) through (20) one obtains 

Im/i = -TImSi , (49) 

where Si is given by Eq. (Il9|) : 

^1 = -^ E Piixi,xi+i)Ln{xi\\^^Ti\xi+i) . (50) 
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With Eqs. ( [T^ ) and (^6|), it is easily shown that 



Ar^Ti = Ar% - io) (>-^)(o|. 



Ai 



Ao 



(51) 



A term in the sum in Eq. ( pOD contributes to the imaginary part of Si only if the argument of 
the logarithm is negative. If we let Y denote the set of all configurations {xi,xi^i) such that 
{xi\\i^Ti\xi+i)<0, then for all {xi,xi^i)(^Y 



Im(Ln(x;|Ai ^Ti|xz+i)) = -ivr . 
Using Eqs. (^9[), (50) and (52) one obtains 



T 



Im/i = -vr— ^Pi(x«,Xi+i) , 



(52) 



(53) 



Y 



where J2y denotes the sum over all configurations {xl,Xl-^-l)€:Y. With Eqs. (10) and ( |5l| ) this can 
be written as 



T 



Im/i = -71"— 



(54) 



Y 



{xi\X]:^To\xi+i) - {xi\0) 



Aq Ai 
Ai Aq 



{0\xi 



+11 



Using the definition of the correlation length in a system of finite width N, ^^^=ln |Ao/Ai|, and 
the fact that ^tv— ^cxd for T^O, one obtains 



Aq Ai 
Ai Ao 



2sinh(e]^i) , 



so that 



Im/i 



(55) 



(56) 



Y 



{xi\X{'To\xi+i) - {xi\0)2smh{Q'){0\xi+i)\ (xz+i|l) 



The quantity in square brackets in Eq. (^) is simply (xi|Aj"^Ti|x;_|_i), and therefore negative 
for every (a;;, x^+i) included in the sum. Since each of the terms in the difference inside the square 
brackets is positive, we must have 



< {xi\Xi^To\xi+i) < (xi|0)2sinh(Gi)(0|x;+i) 
and consequently, since sinh(^^"'^)=^^"'^+0(^^^), we obtain 

vrT 

Im/i = const— + 0(C^3) ^ (iv^^)-i , 



(57) 



(58) 



which yields Eq. (|40|). Only two exceptions to this scenario are possible. Exact cancelation of the 
leading-order terms in the square brackets for all {xi,xi+i)€Y , or multiplicative factors of 0{^^^) 
arising from the matrix elements with the eigenvectors. Both effects are highly improbable and in 
any case would lead to an even smaller |Im/i| than that given by Eq. (|58|). 
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For A^=l, Ti is a 2x2 matrix, so that we can easily calculate Im/i exactly, obtaining 



Im/i 



(59) 



This gives const=— 1/2 in Eq. (58). Under the assumption that the full transfer matrix of the system 
at H=0 and for A^>1 can be approximated by an effective 2x2 transfer matrix with eigenvalues 
that are identical to Aq and Ai of the full transfer matrix |l85| , p^ , 87], one would expect 



Im/i 



--(1 



) 



(60) 



leading to const«— 1/2 for A^>1 as well. In Fig. A.l we plot 2A^|Im/i|/(7rr,^^"'^) as a function of 
for all the temperatures studied. Our expectation is met for low temperatures and small system 
sizes. However, the plotted quantity still depends rather strongly on N and T for higher T and 
N. This is not so surprising, since for increasing T entropy terms also increase, and for increasing 
N the eigenvalue spectrum becomes more complex, both effects affecting the values of the matrix 



elements in the sum in Eq. (56). 

We have observed that the quantity plotted in Fig. A.l shows nice data collapse onto a single 
scaling function when plotted vs. the scaling variable \/N — IT. However, we have not yet found 
a theoretical argument to support this empirical scaling relation. 



B On the applicability of the compact-droplet picture 

For the interpretation of the results of the CTM method in terms of the droplet model to be 
consistent, the critical droplets have to have well defined boundaries. This requires that the size 
of the critical droplets is larger than the single-phase correlation lengths in both the stable and 
metastable phases. A single-phase correlation length, ^gp, can be defined by 



4p' = In 



A, 



a 



A 



(61) 



For the single-phase correlation length corresponding to the stable phase, Aq,=Ao is the dominant 
eigenvalue of the TM, whereas A/j is given by the eigenvalues of the eigenstates that have magne- 
tizations closest to M=(N — 2)/N. In the case of the metastable phase, Aq are the eigenvalues 
corresponding to the metastable branch in the eigenvalue spectrum, whereas A^ are the eigenvalues 
corresponding to the branches with magnetizations closest to M=—{N — 2)/N. 
An estimate for the eigenvalues Aq, can be obtained analytically by using 

N N 

In A, « --Re/a = --([/,- MaH - TReS^) , (62) 

where Ua, Ma, and 5^ are the constrained internal energy, constrained magnetization, and con- 
strained entropy of the eigenstate |a), as defined in Eqs. (p!7[), (p^), and (|l9|). For low temperatures 
we can neglect the entropy term and approximate the internal energy and the magnetization by 
their values at T=0. We therefore insert C/a=— 2J, Up=Ua+{4:J/N), and Mq, and as given 
above into Eqs. (|6ll ) and (^). This yields: 

" 2(2J + H) 
T 

^"'^ ^ 2(2J - H) ■ ^^^^ 
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Here ^st and denote the single-phase correlation length of the stable and the metastable phase, 
respectively. 

Figure B.l shows a comparison between the critical diameter, 2Rc, and the single-phase corre- 
lation lengths, ^st and ^ms) for T/J=1.2. The solid lines were obtained from Eq. (|3^ ) for Rc and 
from Eqs. (S3) and ( |6^ for the single-phase correlation lengths. Also shown are the single-phase 
correlation lengths as obtained from TM calculations for A^=10. For the stable single-phase cor- 
relation length the agreement between the TM results and Eq. ( |63|) is excellent. In the case of 
the metastable single-phase correlation length the TM calculations agree well with Eq. (|6^) up to 
H/J^l. For fields beyond H/J^l, the increasing discrepancy between the metastable single-phase 
correlation length as obtained from the TM and Eq. (^), is caused by entropy effects and the 
deviation of Ua and Mq, from their zero-temperature values. 

An estimate for the magnetic field -ffiviFSP) at which standard continuum nucleation theory is 
expected to become suspect, is given by the intersection of ^ms and 2Rc. An estimate of the temper- 
ature dependence of -ffiviFSP can be obtained by equating from Eq. (|6^ and 2Rc from Eq. (p8|). 
Using the exact zero- field equilibrium surface tension aQ^(T) and the equilibrium magnetization 
meq(r) we get 

gMFSp(T) ^ 4al\T) 

J [2a^^(r) + rmeq(r)] ' ^ ' 

This relation is illustrated in Fig. B.2. As expected, Hmfsp/J=2 at T/J=0. However, since 
Eqs. ( |6^ ) and ( |6^ ) are low-temperature approximations, Eq. (^) must be unreliable as T approaches 
Tc. Certainly, continuum droplet theory breaks down when the diameter of the critical droplet 
becomes smaller than one lattice unit |8^. Using this condition one obtains i^MFSP as shown 
by the dashed line in Fig. B.2. The -f^MFSP obtained from this condition is of the same order of 
magnitude as the one obtained from Eq. (|65|). Also plotted in Fig. B.2 (shown as o) is if™™^(T) 
for A^=10 (from Table |^(c)) for all those temperatures that permitted us to determine it with the 
CTM method. Since i/™™^(T) is the field at which according to Table ||(c) 2i?c~l, we expect it 
to be of the same order of magnitude as ^^mfsp- This is indeed seen to be the case. The results 
summarized in this appendix confirm the applicability of the compact-droplet picture with sharp 
interfaces in the entire field and temperature range covered by our study. 
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Figure 2: (a) The eigenvalue spectrum, —\n\a/{l3JN)^ and (b) the constrained magnetizations 
Mq, as given by Eq. (p^), as functions of the magnetic field H/J for r/J=1.0 {T /TcK,QAA1) and 
N=8. Shown are only those branches that correspond to eigenstates that are symmetric under 
translation and reflection. The segments that contribute to the metastable branch are marked by 
thick lines. 



Figure 3: (a) The real parts, Re/a/J, and (b) the imaginary parts, \lmfa\/ J, of the constrained free 
energies vs. H/J, for T/J=1.0 and N=8. Notice the overall similarity of the stable and metastable 
branches of Kef a/ J to those of the eigenvalue spectrum in Fig. |2|(a). The composite metastable 
branch is marked by thick lines. See detailed discussion in Sec. ^ 



Figure 4: Semi-log plot of the imaginary parts of the constrained free energies, \lmfa\/J, that 
correspond to the metastable branch, shown vs. J/H. Plotted are data for two different system 
sizes, A^=9 and N=10, at T/J=1.0. The thick, straight line was drawn through the two minima 
between J/H=3.0 and J/i?=4.0 (marked by •) for A^=10 only. 



Figure 5: Comparison of the metastable A^|Im/Q,|/J (*) and ^j^^ (o) vs. in a semi-log plot at 
H=0 and T/J=1.0. The dashed straight lines are guides to the eye drawn through the points at 
N=9 and N=10. See details in Sec. O. 



Figure 6: Semi-log plot of the minimum of jlm/aj/J at the weak field 0<ff™™<i7i, plotted vs. 
for T/J=0.4, 0.6, 0.8, 1.0, and 1.2. At these weak fields, llm/d decays exponentially with A^ as is 
evidenced by the straight lines. For each temperature, these lines were simply drawn through the 
points corresponding to the two largest system sizes available. 



Figure 7: Comparison of extrapolations of the surface tension, a{T)/J, as obtained from Eq. 
(shown as x), with the exact equilibrium surface tension, cJeq(T) (solid curve). Also shown (□) are 
estimates obtained by adding a bulk term of the form 2H{N — 1) in the exponential in Eq. (Bot). 



See detailed discussion in Sec. 4.1 



Figure 8: (a) Log-log plot of iff ™/ J vs. A^ for T/J=0.4, 0.6, 0.8, 1.0, and 1.2. The straight hues 
show that i/f™ is proportional to N~°'^'^\ (b) The exponent a(T) as a function of T, as obtained 
from linear least-squares fits to our data. 



Figure 9: Shown is a plot of the two-point finite-difference estimate for /3H(T)/J, 
-A(ln|Im/„|)/A(J/|iJ|), as a function of H^s/J (compare Eqs. (|^ and (H)) for four differ- 
ent temperatures. The symbols represent the results obtained with the CTM method, except at 
zero field, where we plotted /3Heq(T)/J obtained from the exact equilibrium quantities Seq(T) and 
Am=2meq(r) as described in Sec. [4.2| . The solid curves represent three-parameter linear least- 
square fits. The dotted straight lines have a slope of unity, corresponding to a value of 6=1, and 
they intercept the vertical axis at the exact /3Heq(T)/J. 
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Figure 10: The quantity H(T) divided by its equilibrium value Heq(T), calculated from a three- 



parameter linear least-squares fit to Eq. (47). The O correspond to results from fits including all 
the minima of the metastable |Im/Q,| in the interval i7™°</f<2J, whereas the x with error bars 
were obtained by excluding to reduce uncontrollable finite-size effects. (For clarity, no error 

bars are shown for the O.) The system size was A^=10, except at T/J=0.4, where we used A^=9. 



See Sec. 4.2 for a detailed discussion. 



Figure 11: The prefactor exponent b vs. T/ J. The values shown were obtained from a two-parameter 
linear least-squares fit to Eq. (^) by setting H(r) to its equilibrium value. See Sec. L2 for details. 



Figure 12: The quantity H(T)/J^ as obtained from a two-parameter linear least-squares fit to 
Eq. (p7|), setting 6=1. The solid line corresponds to the equilibrium value Heq(r). The droplet 
shape interpolates between a square at T=0, given by (2(TQ)/meq, and a circle at T=Tc, given by 
(7r(To)/(2meq), both of which are shown for the whole temperature range as dashed curves. The 
metastable CTM estimates, shown as • with error bars, follow the equilibrium curve closely. The 



dashed vertical line marks the critical temperature. See Sec. 4.2 for details 



Figure 13: Semi-log plot of the prefactor B{T) vs. inverse temperature J/T. The CTM results are 
indicated by x with error bars and were obtained from a two-parameter linear least-squares fit to 
Eq. ( p9|) by setting H(T) to its equilibrium value and 6=1. The dashed line is a weighted linear 
least-squares fit of the values oi B{T) to the form B{T)=kexp{-AF/T), and gives AF/J=1.78(5) 
and /c=2.12(7). Also shown are Monte Carlo results for B{T) (O with error bars) as obtained by 
continuum droplet theory from direct measurements of the lifetime of the metastable phase ||5^, ^ . 
The dot-dashed line is a weighted linear least-squares fit which gives AF/J=1.25(5) and /c=3.54(8). 
The o represent B{T) as calculated from Becker-Doring droplet theory by Harris ||37[| . See Sec. [4.2| 
for details. 



Figure .1: The quantity 2N\lmfi\/{7rTQ'^) at H=0 function of N is shown for all the tem- 
peratures studied. At T/J<0.8 numerical underfiow prevented us from obtaining reliable results 
for the larger N. An exact calculation of 2A^|Im/i|/(7rT^]^"^) at H=0 for A'^=l leads to a value of 
unity for all temperatures (shown as •). See Appendix]^ for more details. 



Figure .1: Comparison of analytical estimates of the diameter of the critical droplet, 2i?c, with 
the estimated single-phase correlation lengths, .^st for the stable phase, and ^ms for the metastable 
phase at T/J=1.2. Here 2Rc was calculated from Eq. (^), and ^st and ^ms were obtained from 
Eqs. (]6^ ) and (|6^). The transfer-matrix results for and ^ms foi' A^=10 are shown as O. 



Figure .2: Estimates of the crossover field Hmfsp/J obtained from Eq. (65) (solid line) and from the 
condition 2Rc=l (dashed line). Also plotted (o) is HY^™-^{T)/J for A^=10, the field corresponding 
to the minimum of the metastable llm/d in the field interval 1<H/ J<2. 
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